/*! \file spai.h
 *  \brief Sparse Approximate Inverse (SPAI) preconditioner.
 */

#pragma once

#include <cusp/linear_operator.h>
#include <cusp/csr_matrix.h>


namespace cusp
{
namespace precond
{


/*! \addtogroup preconditioners Preconditioners
 *  \ingroup preconditioners
 *  \{
 */

/*! \p Sparse Approximate Inverse
 */
template <typename ValueType, typename MemorySpace>
class spai : public linear_operator<ValueType, MemorySpace>
{
	typedef linear_operator<ValueType, MemorySpace> Parent;

public:

	typedef typename cusp::csr_matrix<int, ValueType, MemorySpace> matrix_type;

	matrix_type matrix;

public:
	/*! \brief Construct a spai preconditioner.
	 * Construct an approximate inverse of \p A, using \p A as sparsity pattern.
	 *
	 * A must be a square matrix, preferably in cusp::csr_matrix format.
	 *
	 * \param A matrix to invert
	 */
	template<typename MatrixType>
	spai(const MatrixType& A);

	/*! \brief Construct a spai preconditioner.
	 * Construct an approximate inverse of \p A, using \p P as sparsity pattern.
	 *
	 * A must be a square matrix, preferably in cusp::csr_matrix format.
	 *
	 * \param A matrix to invert
	 * \param P sparsity pattern for the preconditioning matrix
	 */
	template<typename MatrixType, typename MatrixType2>
	spai(const MatrixType& A, const MatrixType2& P);

	/*! \brief Apply the preconditioner to vector \p x and store the result in \p y
	*
	* \param x input vector
	* \param y ouput vector
	* \tparam VectorType1 vector
	* \tparam VectorType2 vector
	*/
	template <typename VectorType1, typename VectorType2>
	void operator()(const VectorType1& x, VectorType2& y) const;

private:
	template<typename MatrixType, typename MatrixType2>
	void make_spai(const MatrixType& A, const MatrixType2& P, cusp::host_memory, cusp::known_format, cusp::known_format);

	template<typename MatrixType, typename MatrixType2>
	void make_spai(const MatrixType& A, const MatrixType2& P, cusp::device_memory, cusp::csr_format, cusp::csr_format);

	template<typename MatrixType, typename MatrixType2>
	void make_spai(const MatrixType& A, const MatrixType2& P, cusp::device_memory, cusp::known_format, cusp::known_format);

};

/*! \}
 */

} // end namespace precond
} // end namespace cusp

#include <cusp/precond/detail/spai.inl>

